##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/DriverChoose/GeneVariantType"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_info <- opt$sample_info
out_path <- opt$out_path
ccf_file <- opt$ccf_file

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_info <- data.frame(fread(sample_info))
dat_ccf <- data.frame(fread(ccf_file))

###########################################################################################

Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
gene_list <- c("TP53")

###########################################################################################
## 加载分类
dat_ccf <- merge( dat_ccf , unique(dat_info[,c("ID" , "Type")]) , by = "ID" )
dat_ccf <- subset( dat_ccf , Type != "IM + IGC + DGC" )

dat_ccf2 <- dat_ccf
dat_ccf2$Type <- "All"
dat_ccf <- rbind( dat_ccf , dat_ccf2 )

###########################################################################################
## 判断基因多少比例出现在Trunk、Pre_Private、Inv_Private
## 计算p值
trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

result_plot <- c()
for(typeN in unique(dat_ccf$Type)){
    print(typeN)
    for( geneN in gene_list ){
        print(geneN)
        tmp <- subset( dat_ccf , Variant_Classification %in% Variant_Type & Hugo_Symbol == geneN & Class != "IM" & Type == typeN )
        tmp$Variant_Classification <- ifelse( tmp$Variant_Classification %in% c("Nonsense_Mutation" , "Frame_Shift_Del" , "Frame_Shift_Ins" , "Splice_Site") , "LOF" , tmp$Variant_Classification )
        tmp$Variant_Classification <- ifelse( tmp$Variant_Classification %in% c("In_Frame_Del" , "In_Frame_Ins") , "In_Frame" , tmp$Variant_Classification )
        tmp$Variant_Classification <- ifelse( tmp$Variant_Classification %in% c("Missense_Mutation") , "Missense" , tmp$Variant_Classification )

        res_tmp <- c()
        for( sample in unique(tmp$ID) ) {
            tmp_use <- subset( tmp , ID == sample )
            tmp_use <- unique(tmp_use[,c("ID" , "Share" , "Variant_Classification")])

            ## 一个人若发生多个突变，算共享的的
            if( length(which(tmp_use$Share==TRUE)) > 0 ){
                tmp_use <- subset( tmp_use , Share == TRUE )
            }

            if( nrow(tmp_use) > 1 ){
                tmp_use <- tmp_use[1,]
            }

            res_tmp <- rbind(res_tmp , tmp_use)
        }

        res_tmp2 <- data.frame(table(res_tmp$Variant_Classification , res_tmp$Share))
        p <- fisher.test(matrix(res_tmp2$Freq , ncol = 2))$p.value
        
        if( p < 0.001 ){
            p_text <- trans(p)
        }else{
            p_text <- paste0( "P == " , round(as.numeric(p) , 2) ) 
        }

        res_tmp2$p <- ""
        res_tmp2$p[1] <- p
        res_tmp2$p_text <- ""
        res_tmp2$p_text[1] <- p_text
        res_tmp2$Gene_Symbol <- geneN
        res_tmp2$Type <- typeN
        res_tmp2 <- res_tmp2 %>%
            group_by(Var2) %>%
            mutate(count_all=sum(Freq),
               ratio=Freq/count_all)

        result_plot <- rbind( result_plot , res_tmp2 ) 
    }
}

out_name <- paste0(out_path , "/VariantType.TP53.tsv")
write.table( result_plot , out_name , row.names = F , sep = "\t" , quote = F )

###########################################################################################
## 基因分布堆叠图
result <- data.frame(result_plot)
result$ratio[is.na(result$ratio)] <- 0
result$value_text <- round(result$ratio , 2) * 100
print(result)

result$Var1 <- factor( result$Var1 , levels = c("Missense" , "In_Frame" , "LOF" ) )
result$TrunkUse <- ifelse( result$Var2 == "TRUE" , "Trunk" , "Private" )
result$value_text <- ifelse( result$ratio == 0 , "" , result$value_text )
result$TrunkUse <- paste0(result$TrunkUse , "(" , result$count_all , ")")
result$TrunkUse <- factor( result$TrunkUse , levels = unique(result$TrunkUse )[order(unique(result$TrunkUse ) , decreasing=T )] )

## 提取TP53
plot <- ggplot( data = subset(result , Gene_Symbol == "TP53") , aes( x = TrunkUse , y = ratio , fill = Var1 ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Proportion (%)")+
    facet_grid(.~Type,space='free_x',scales='free_x') +
    theme(panel.grid = element_blank())+
    scale_fill_npg() +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1.05 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='right',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 8,color="black",face='bold'),
                axis.text.y = element_text(size = 7,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.text.x = element_text(size = 8,color="black",face='bold') ,
                axis.line = element_line(size = 0.5))

out_name <- paste0(out_path , "/VariantType.TP53.pdf")
ggsave( out_name , plot , width = 6 , height = 4 )
